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We propose a simple nonlocal energy functional that is suitable for the continuum field charac- 
terization of nonperiodic and localized textures. The phenomeno logical functional is based on the 
pairwise direction-dependent interaction of field gradients that are separated by a fixed distance. In 
an appendix, we describe the numerical minimization of our functional. On that basis, we investigate 
the kinetic evolution of thread-like stripe patterns that are created by the functional when we start 
from an initially disordered state. At later stages, we find a coarse- graining that shows the same 
scaling behavior as was obtained for the Cahn-Hilliard equation. In fact, the Cahn-Hilliard model 
is contained in our characterization as a limiting case. A slight modification of our model omits this 
coarse-graining and leads to nonperiodic stripe phases. For the latter case, we investigate the tem- 
poral evolution of the defects (end points) of the thread-like stripes. In view of possible applications 
of this functional, we discuss a possible characterization of polymeric systems and vesicles. The 
statistics of the growth of the thread-like structures is compared to the case of step-growth poly- 
merization reactions. Furthermore, we demonstrate that the functional may be applied for the study 
of vesicles in a continuum field description. Basic features, such as the tendency of tank-treading in 
simple shear fiows and parachute folding in pipe fiows, are reproduced. 

PACS numbers: 81.16.Rf,83.80.Qr,87.16.D-,81.16.Dn 



I. INTRODUCTION 

For several decades, concentration and density field de- 
scriptions have been used to model the behavior of pe- 
riodic micro- and mesoscopic systems. A prominent ex- 
ample is given by the Ohta-Kawasaki characterization of 
mesophases as they are observed in block copolymers [T] . 
This description was derived by averaging over the local 
positions of the monomers that form the different blocks 
of the polymers. The procedure led to an effective contin- 
uum characterization of the observed periodic structures 
by the famous Ohta-Kawasaki energy functional. In the 
case of incompressible diblock copolymers, the order pa- 
rameter is given by the concentration of one of the two 
monomer types that form the different blocks. This is 
why we call this sort of characterization a concentration 
field description in the following. 

Since then, a lot of studies have demonstrated the ef- 
fectiveness of the approach. For instance, the kinetic 
evolution of the periodic mesophase textures from the 
disordered (spatially homogeneous) state has been inves- 
tigated 0. Also the kinetic evolution from one texture 
to the other after temperature quenches has been ana- 
lyzed [3l H] . Elastic moduli characterizing deformations 
of the textures were determined [5^ and rheological 
properties of such structures were studied [THTT], 

In a recent model of similar spirit, an energy functional 
has been postulated to characterize the behavior of crys- 
talline structures through a density field description. We 
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refer to the phase field crystal model suggested by El- 
der and Grant [T2[ [13] . Periodicity of the density field 
is put into the model by energetically stabilizing at least 
one nonzero wave vector in the ground state. The mag- 
nitude of the subsequent periodic density field is inter- 
preted as the averaged probability to find a constituent 
of the crystal. Short timescales of motion are regarded 
as averaged over. Because of this, the approach is very 
effective in modeling processes in crystals on timescales 
that are long from a molecular dynamics point of view. 
We refer to such a model as a density field description 
because the order parameter is related to a mass density 
rather than a concentration field. 

Kinetic evolution of periodic patterns from a disor- 
dered state and corresponding phase diagrams have been 
analyzed in the framework of the phase field crystal 
model [13 . Furthermore, the approach has been used 
to describe the motion of defects in crystalline structures 
[14, 15 , including avalanche events [16 . The phase field 
crystal model has been shown to be a starting point for 
multiscale simulations of multidomain structures [171 [18]. 
In addition, the suitability for the characterization of epi- 
taxial growth and crack propagation has been pointed 
out [m [13]. It has been demonstrated that averaging 
the trajectories obtained from molecular dynamics simu- 
lations indeed leads to density fields that can be captured 
by the phase field crystal model [19 . The connection to 
classical density functional theory is of current interest 
[201 [2l. 

Both of the two approaches outlined above can be 
interpreted as describing a sort of phase separation on 
micro- or mesoscopic length scales. This means that 
macroscopic coarse-graining of the observed structures, 
which sets in, for example, in the case of the Cahn- 
Hilliard model [22 , is suppressed. The structures that 
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are obtained are periodic in space. Or, at least, if we 
think of defect-rich and/or multi-domain structures, the 
energetically favored ground state is periodic in space. It 
is our goal in this paper to describe a sort of microphase 
separation that does not feature this tendency and is 
nonperiodic in space. Whereas the phase field crystal 
model has been successful to characterize processes in 
crystalline materials, our approach may be interesting to 
model amorphous materials such as polymers. 

In the next section we phenomenologically introduce 
such a nonlocal energy functional and motivate the con- 
cept behind it. We explain how we numerically solve the 
corresponding kinetic equation in the appendix. After 



that, in section III, we give an overview on the textures 
that are obtained from minimizing the energy functional. 
In particular, we focus on the time evolution of the struc- 
tures as they develop from the disordered state. At this 
stage, coarse-graining to macrophase separation still sets 
in. We therefore present a modification to our model in 
section [TVl which leads to stabilized stripe phases that do 
not phase separate macroscopically. Some possible appli- 
cations, predominantly in the field of soft matter systems, 
are pointed out in section [V| before we conclude. 



II. SIMPLE NONLOCAL ENERGY 
FUNCTIONAL 

As pointed out above, we now introduce a simple 
energy functional that produces localized nonperiodic 
structures. The functional is introduced in a phenomeno- 
logical way. However, we will shortly explain the concept 
behind the chosen functional form using the illustration 
in Fig. [l] After that, we will point out the sort of tex- 
tures obtained. We will restrict ourselves to two spatial 
dimensions in this paper. 

Our goal is to use a single scalar order parameter field 
ij{r) to describe localized nonperiodic objects, as for ex- 
ample the ones shown in Fig. [ija). tp{r) denotes a con- 
centration or density field, similar to the continuum field 
description on block copolymers pLj and the phase field 
crystal model [13]. The spheres and stripes shown in 
Fig.jlja) feature a characteristic size d that we externally 
impose. It corresponds to the diameter and thickness of 
the spheres and stripes, respectively. In this example, 
the distances between the objects do not follow a peri- 
odic rule. More details are given below. 

Without loss of generality, we identify the black and 
white regions in Fig. [TJa) with values ip = —1 and 

= +1, respectively. Fig. [TJc) shows the density or 
concentration profile along a horizontal cut through one 
of the stripes. 

We now consider a single isolated up-gradient, as 
shown, for example, in the left half of Fig. [TJc). It would 
be energetically unstable in periodic approaches like the 
ones introduced above [H [13] because it is not embedded 
in an overall periodic structure. Our idea is to balance 
each such local up-gradient by a down-gradient that is 
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FIG. 1: (Color online), (a) Textures observed from a numer- 
ical iteration of Eqs. ([T]) and ^ forward in time, starting 
from a homogeneous disordered state. Brightness indicates 
the value of ip, where black and white correspond to values 
around —1 and +1, respectively. The actual calculation do- 
main was of size 1000 x 20. Only a domain of size 150 x 20 
is shown. The boundary conditions are periodic, (b) Density 
profiles '02,j=i5 along a horizontal line at height j = 15 in 
Fig. (a) (+ symbols) and ip^^j averaged over the j-direction 
of Fig. (a) (x symbols), (c) Zoom into the area in (b) sur- 
rounded by the dashed rectangle. Gradients Vip at the lattice 
sites pairwisely sum up to zero, as indicated for one pair. The 
distance between the points of each such pair is A = 14 lattice 
distances. (Technical details: = 5, a = 3, d = 15, L = 1, 
ip = —0.4, lattice distance = 1, 800000 time steps of step 
size dt = 0.0008 yielding a time t = 640). 



located a distance d away in the direction the gradient 
is pointing to. This direction is given by the unit vector 
V?/^/|| VV^||. In this way, we form localized pairs of gra- 
dients that energetically balance each other, having the 
same magnitude but opposite sign and being separated 
by an imposed distance d A in Fig. [TJc), see below). 
Translating this into a formula, the energy functional 
that we will investigate is based on the energy density 



F(r) 



1 



(1) 



Here, 'i^ is a measure for the inverse temperature. The 
part in square brackets is weighted by a second constant 
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parameter a and contains the impact of the gradient field 
VV^. Here, the first term refers to at position r and 
is local. The second term, however, refers to the value 
of S/ip at a position of distance d away from r in the di- 
rection that points to. This latter term is therefore 
a nonlocal contribution. Expression ([T]) satisfies the nec- 
essary symmetry requirements. To guarantee thermody- 
namic stability, we must set > and a > 0. We obtain 
the energy functional J^{'0(r)} via T = J F{t) d?r. 

The kinetic equation for the time evolution of the con- 
served field tp{v) follows from minimizing the energy J^, 

^=LA^. (2) 

dt Si; ^ ^ 

Mostly in this paper, the parameter L will be set equal 
to one, L = 1, which can always be achieved by rescaling 
the time. On the right-hand side of this equation, the 
Laplacian refiects the fact that ip is a conserved quan- 
tity. The nonlocal contribution of Eq. (fTl) is implicitly 
contained in the functional derivative dT/dip. 

We must remark at this point that, through the form of 
Eq. ([2|, we limit our investigations to the regime of relax- 
ation dynamics. Processes on time scales that are shorter 
than the diffusive one have to be considered as being av- 
eraged over. Such processes cannot be resolved explicitly 
by this approach. As a benefit from this coarse-grained 
procedure, however, longer time scales are accessible by 
numerical calculations. The same approach was used in 
the introduced studies on block copolymer systems [Ji- 
ll] and the phase field crystal model [IlHIl Ull UMII], 
or, in a similar spirit, also in classical dynamical density 
functional theory [2ni23H25]. 

Due to the nonlocal contribution, we could not per- 
form the functional derivative on the right-hand side of 
Eq. ([2| analytically. We present our way of numerical 
evaluation of the functional derivative in the appendix 
and the supplemental material [26]. The central point in 
our discretized numerical evaluation is to track the con- 
nections between the points r and r + d^^. Due to the 
discrete nature of the numerical results, we will replace 
the continuous position vector r by indices i and j. 

In the current form, Eq. ([T]) contains three parameters 

a, and the distance d. It follows directly that a can 
be scaled out. For practical reasons, we set a = 3 in all 
our numerical considerations. We further note that in 
the limit (i ^ 0, our expression for F leads to the famous 
Cahn-Hilliard equation [22 . 

A first illustrative example created from our nonlocal 
energy functional is actually given by Fi g. \\\ It shows 
the result of an iteration of Eqs. ([T]) and ([2| forward in 
time on a rectangular domain of large aspect ratio and 
periodic boundary conditions. The axis labels i and j 
index the grid sites of the square calculation grid. We 
initialized the field ip as ipij = + 77^^, with ?/5 = —0.4 
and Tjij a random field of zero mean, Gaussian distribu- 
tion, and small amplitude. Sphere and stripe domains 
are obtained. Their diameter follows from the value of 
the parameter which we set to d = 15. 



We display in Fig. [TJb) the density profile tpij=i5 on 
a horizontal line at j = 15, and a vertically averaged 
profile ip^j = Xlj=o ^^j/^-'-- "^^^ values of the two pro- 
file functions are equal for the stripes and differ for the 
spheres. 

Fig. [TJc) is a zoom into the area of Fig. [l{b) that IS 
marked by the dashed rectangle and depicts the concept 
behind Eq. ([T]). The first part with the coefficient i9 en- 
forces concentrations of ip = ^1 for > 0. All other val- 
ues increase this contribution to the energy of the system. 
The second part with the coefficient a leads to an ener- 
getic penalty for gradients V?/^, unless there is a gradient 
of opposite sign, but same magnitude, located at a dis- 
tance d away in the direction that \/tp points to. Indeed, 
this pairwise annihilation of gradients occurs as indicated 
for one pair in Fig. [ijc). (The distance is A = 14 and 
not d = 15 because of an integer cast of the term 
in our numerical calculation to index the discrete lattice 
sites). In this way, no energetic penalty arises from the 
gradient terms in Eq. ([T]). 

In summary, we do not confine ourselves to local ex- 
pressions in the gradient fields. Such local expressions 
have the tendency to enforce periodicity throughout the 
whole space. This becomes clear from the Fourier trans- 
forms of the corresponding energy densities. Usually, a 
certain band of wave numbers is preferred homogeneously 
over the whole space. Our goal is to generate objects that 
are localized in space and not embedded in an overall 
tendency of periodicity. This is different from multigrain 
structures or defect-rich textures since in the latter cases 
a periodic order is still energetically favored. In other 
words, we wish to enforce phase separation on a meso- 
scopic length scale that leads to amorphous textures as 
observed in non-crystallized polymeric systems, for ex- 
ample, or localized objects such as vesicles ffoating in a 
ffuid environment, using one scalar order parameter field. 



III. TIME EVOLUTION OF THE TEXTURES 

In the following, we describe the kinetic evolution of 
the patterns that are generated by Eqs. ([T]) and Since 
we compare to results obtained for the Cahn-Hilliard 
equation, we refer to the order parameter as a concen- 
tration field. We will start our numerical considerations 
in all cases from a spatially homogeneous concentration 
that is modulated only by a random field of zero mean, 
Gaussian distribution, and small amplitude. The numer- 
ical calculations were performed in the way pointed out 
in the appendix and the supplemental material [26 . 

We show in Fig. [2] a time series of phase separation 
for a mean concentration tp = 0. Furthermore, we chose 
I? = 5, a = 3, (i = 5, and a grid size of 512 x 512 lat- 
tice sites. The brighter the color in the pictures, the 
higher the concentration (for each picture the concentra- 
tion values were normalized by their maximum value). 
Fig. [3] represents the analogous time series for a mean 
concentration ip = —0.4. 




FIG. 2: Time series of the kinetic evolution described by Eqs. ([l]) and ([2|, with a mean concentration -0 = 0. The gray scale 
represents the concentration values, where white and black correspond to the maximum and minimum concentration at each 
time, respectively. Parameter values were set to = 5, a = 3, d = 5, the grid size to 512 x 512 sites of distance dx = 1. Time 
steps were of size dt — 0.0008, the times of taking the pictures are indicated in their bottom right corners. 



First, thread-like stripe structures form. Shorter 
threads connect to longer threads. The reason for these 
steps of forming longer objects is the higher energy den- 
sity at the ends of the threads: there is no concentration 
gradient within the threads that could balance the lon- 
gitudinal concentration gradient at the end points. This 
increases the local energy density (compare Eq. ([T]) and 
Fig. [T]). The upper part of Fig. |4] highlights the spatial 
distribution of the energy density at an early stage of pat- 
tern evolution for the case = 0. Ends of the threads 
clearly stand out. We further investigate this initial pro- 
cess of forming and connecting to longer threads in the 
next section. 

After that, another process of coarse- graining sets in, 
which leads to macroscopic phase separation. Blobs start 
to form at defects in the thread structure and then pref- 
erentially grow and coarse-grain along the thread con- 
nections. The smaller the mean concentration the less 
connected are the blobs by threads during the process 
of coarse-graining. As we will show below, this coarse- 



graining process is similar to the one observed for the 
Cahn-Hilliard model. 

Energetically, we can explain this coarse-graining to a 
macroscopic phase separation in the following way. The 
stripe patterns as they evolve from Eqs. ([T]) and ([2| corre- 
spond to a metastable state. This is true even if the high- 
energetic end points of the threads are excluded from our 
considerations. Due to pairwise annihilation of the gradi- 
ents at the boundaries of the stripes, the term with coeffi- 
cient a can be minimized in Eq. ([T]). However, the values 
at the boundaries of the stripes deviate from = ±1 
(compare Fig. [T]), which increases the '^^-contribution. 

For the blobs, the situation is just the other way 
around. The gradients at the surfaces cannot be anni- 
hilated, so the a-term in Eq. ([T]) increases in magnitude. 
This leads to the high energy density at the surfaces de- 
picted in the lower part of Fig. [4] Yet, this scenario is 
energetically beneficial. In most of the interior of the 
blobs we find ^ 1, which implies F ^ 0. Altogether, 
this decreases the total energy of the system due to the 
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FIG. 3: Time series of the kinetic evolution described by Eqs. ([T]) and ([2]), with a mean concentration -0 = —0.4. The other 
settings were the same as described in the caption of Fig. |2] 



high number ratio of interior to surface points. 

To characterize the coarse-graining process more quan- 
titatively, we determined the structure function 



structure function 



S(k,() 



-ikr 



(3) 

Here, N is the total number of lattice sites. The sums run 
over all lattice sites (r and correspond to index pairs 
(z, j) and m), respectively, that index the lattice sites). 
Likewise, the vector k refers to discrete lattice sites of the 
corresponding discrete reciprocal lattice. (...) denotes 
an ensemble average. In our case, we average over three 
independent numerical runs that only differ in the initial 
conditions. More precisely, the three numerical calcula- 
tions start from different spatial realizations of the nar- 
row Gaussian distribution of concentrations around the 
spatially homogeneous state. 

From 6'(k, t), we calculate the circularly averaged 



k' < k-\-dk 

S{k,t) = - Yl 



k'>k 



(4) 



# stands for the number of lattice sites that satisfy 
the requirement k < < k -\- dk, where we set dk = 2. 
Despite the anisotropy in our numerical results, the circu- 
larly averaged S{k,t) is suitable for what we demonstrate 
below. 

In Fig. |5j we plot the structure function S{k^ t) that 
corresponds to the case of ^5 = and d = b (compare the 
time series in Fig. [2|. Fig.jsja) shows S{k^t) at the early 
times, when the stripe texture leads to a peak around 
/c ^ 65 (we do not display the higher order peaks). This 
peak indicates a periodic length of about 512/65 ^ 8 
lattice distances dx of the thread-like stripe structures. 
Comparing with the first row of pictures in Fig. [2j we 
confirm this length scale: the threads have a thickness 
of d = 5, and, where they are compactly packed, are 
separated by roughly a distance d/2^ which gives d + 
d/2 ^ 8. As time proceeds, this peak becomes smaller 
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FIG. 4: (Color online). Energy density F calculated according 
to Eq. ([T]) for two of the snapshots presented in Fig. [5] the 
ones taken at time t = 20 and at time t — 8000. The reddish 
color marks regions of high energy density. Clearly, the areas 
around defects and at the surface of the blobs stand out. The 
pictures show averages over 100 iterations after the given time 
to smooth the fluctuations that arise from the discrete nature 
of the numerical implementation. Only a domain of 170 x 128 
lattice sites is shown. 



until it has vanished at t ~ 800. 

As the peak at /c ~ 65 decreases, another peak at 
smaller /c-values evolves. This peak corresponds to the 
onset of the second stage of coarse-graining. Fig. [5|b) il- 
lustrates the development of Sik^t) at times longer than 
t = 800. The peak shifts to smaller /c-values, indicat- 
ing a growth of the characteristic size of the observed 
structures. In particular, this is reflected by the last row 
of pictures in Fig. |2] and indicates a macroscopic phase 
separation. 

To discuss the latter statement in more detail, we 
demonstrate that the same scaling behavior is found as 
for the later stages of coarse-graining in the Cahn-Hilliard 
model ^22| ,27| . We follow Ref. ^27^ and define a normal- 
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FIG. 5: (Color online). Circularly averaged structure function 
S{k,t) for a texture evolved from Eqs. |l]) and ^ for ijj = 0. 
The curves correspond to times t as indicated on the upper 
right of each plot. See Fig. [2] for a real space representation 
and the parameter values, (a) Temporal evolution of S{k, t) 
during the early stage. The decreasing first order peak k ^ 
65 represents the thread-like stripe pattern, the increasing 
peak at lower wave numbers develops due to the macroscopic 
coarse-graining. Only data S{k^ t) < 100 are shown, (b) 
Evolution of S{k,t) at the later stage of coarse- graining. As 
the blobs grow in size, the peak shifts to lower /c- values. The 
dashed lines are guides to the eye only, (c) Data collapse for 
the later stage of coarse-graining (see text and Ref. [27|). 
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FIG. 6: Time series of the kinetic evolution described by 
Eqs. ([l]) and ([2]), with a mean concentration = —0.7. The 
other settings were again the same as described in the caption 
of Fig. [2] However, a nucleation core of 5 x 5 lattice points and 
-0 = 0.5 was initially provided in the center of the calculation 
grid. 



ized structure function 



as well as a characteristic wave number 



hit) 



(5) 



(6) 



When we plot ki{t)'^s{k,t) as a function of k/ki{t), the 
data for the later times collapse as displayed by Fig.jsjc). 

After the analysis of this scaling behavior, we further 
increased the magnitude of the mean concentration 
For values > 0.6, we found that the patterns did 
not form spontaneously. This can also be understood 
from a comparison with the Cahn-Hilliard equation. In 
the Cahn-Hilliard model, an initially homogeneous state 
only becomes linearly unstable for mean concentrations 
IV^I < l/>/3 ^ 0.577. We therefore provide a "condensa- 
tion nucleus" at these concentrations > 0.6. Fig. |6] 
shows an example, where we set t/j = —0.7. The other 
parameters were chosen the same as in Figs. [2] and jS] 
However, at t = 0, a spot of 5 x 5 lattice points in the 
center of the figure was set to tjj = 0.5. Indeed, we then 
observe a process of growth after providing the initiat- 
ing nucleation core. Coarse-graining starts at the most 
unstable points in the threads of Fig. |6] 

Finally, we display in more detail an example of how 
defects can form, which then initiates the coarse-graining. 
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FIG. 7: Snapshots of 57 x 57 lattice sites from a numerical 
integration of Eqs. ([T]) and ([2]), with a mean concentration 
'0 — —0.4. We set the parameter values to -i^ = 5, a = 3, 
d = 15, the total grid size was 512 x 512 sites of distance 
dx — 1^ and time steps were of size dt — 0.0008. The times 
t that the different pictures correspond to are (a) 5760, (b) 
6400, (c) 6480, (d) 6760, (e) 7200, and (f) 8000. 



The pictures in Fig. [T] show snapshots from a numeri- 
cal calculation of a mean concentration ip = —0.4 and 
a stripe thickness d = 15. The series starts from two 
threads (a), the energy density at the ends of which is 
increased as noted above. First, the "reactive" end of the 
shorter thread couples into the side of the other thread 
(b). The resulting structure (c) is frustrated since the re- 
quirement of pairwise annihilation of the gradients can- 
not be satisfied around the connection point. This higher 
energy density around the crosslinking point initiates the 
coarse-graining (d). A growing blob incorporates the re- 
maining thread texture (e). Finally, we end up with a 
spherical object (f) that looks reminiscent of a small vesi- 
cle. 

We shortly remark at the end of this section, that the 
preferentially diagonal, horizontal, and vertical orienta- 
tion of the threads is due to our discretization scheme 
of the gradient field as described in the appendix. It re- 
flects the underlying square lattice that acts similar to 
an external orienting field. Varying the value of the pa- 
rameter (i, we checked that the coarse-graining is not due 
to the interplay with the discrete numerical grid. Still, 
the amorphous and localized character of the structures 
is evident, in particular for lower mean concentrations ip 
(Figs. [3]and[6|. We checked the numerical scheme by 
varying the size of the time steps dt^ the distance of the 
lattice sites dx^ the characteristic length d (also to non- 
integer values), as well as the discretization scheme of the 
Laplacian in Eq. ([2|. 

IV. STABILIZED STRIPE PHASES 

To make our model more attractive for the descrip- 
tion of amorphous and localized structures, we need 
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to stabilize the stripe state against macroscopic coarse- 
graining. This becomes clear when we consider surfac- 
tant molecules as an example. The white regions in 
Figs. [2] and |3] are identified with the hydrophobic compo- 
nents, the black regions with the hydrophilic components. 
Then macroscopic phase separation is not allowed. The 
surfactant molecules cannot simply chemically split their 
hydrophilic from their hydrophobic parts. At the same 
time, each gradient in ip is associated with an accumu- 
lation of the hydrophilic component on the up-side and 
the hydrophobic component on its down-side. 

In the following minimal approach, we only consider 
the accumulation on the down-side of each gradient. We 
do this by adding an additional contribution Ff^ to the 
energy density in Eq. ([T]), 



Fp{v) = p i'l.v^ + 1 



TWIT 



(7) 



Due to this contribution, each gradient VV^ favors a value 
of ~ — 1 in its downhill direction. More precisely, this 
term leads to an energetic penalty at the position r, if 
the value of tjj deviates from = — 1 at a distance b 
away in the direction ~w^|r- The interplay with the 
a-contribution in Eq. ([T]) energetically disfavors macro- 
scopic coarse-graining. 

We can understand the latter statement by looking at 
the fine structure of the blob boundaries in Figs. [2| [3j |6| 
and[7| The a-term in Eq. ([T]) enforces the formation of 
concentration shells around the blob, which corresponds 
to the scheme depicted by Fig. [l] It leads to concentra- 
tion gradients at the inner shell boundary pointing into 
the inward direction. However, the inside concentration 
of the blob is i/j ^ +1, which leads to a high energetic 
penalty due to the /^-contribution in Eq. There- 
fore, macroscopic coarse-graining is energetically inhib- 
ited. Details of the numerical implementation of the (3- 
term are given at the end of the appendix and in the 
supplemental material [26 . 

Figs. [S] and [9] illustrate the stability of the stripe pat- 
terns when Eq. Q is included. The mean concentrations 
were given by = and ip = —0.4, respectively. We 
set the other parameter values toi? = 5, a = 3, (i = 5, 
/3 = 1, and b = 3. Again, the grid size was 512 x 512 
lattice points, and the length of the numerical iteration 
was the same as the one leading to Figs. [2]and[3j 

A macroscopic phase separation via the blob forma- 
tion does not occur, despite the presence of defects and 
crosslinking points. The only sort of coarse-graining is 
the formation of longer threads from shorter ones. Two 
different realizations of this process are observed. On 
the one hand, mainly short threads can diffuse into other 
threads. This happens predominantly at the early stage 
of coarse-graining. It requires a nearby thread that can 
grow on the cost of the dissolving thread. On the other 
hand, the "reactive" ends of two different threads can 
combine to form one resulting long thread. If such a 
"reactive end" attacks the inner part of another thread 
similarly to the process depicted in Fig.[7]^b), crosslinking 




FIG. 8: Stability of the stripe state during a kinetic evolution 
described by Eqs. ([T]) and Q, supplemented by Eq. ([7|, for a 
mean concentration -0 = 0. In contrast to Fig. [2] no macro- 
scopic coarse-graining occurs. Parameter values were set to 
1^ = 5, a = 3, d = 5, /3 = 1, and 6 = 3. The grid size was 
512 X 512 lattice points of distance dx — 1^ and time steps 
were of length dt = 0.0008. 




FIG. 9: Stability of the stripe state during a kinetic evolution 
described by Eqs. ([T]) and ([2]), supplemented by Eq. ([7|, for 
a mean concentration '0 = —0.4. No macroscopic coarse- 
graining occurs, in contrast to Fig. |3] Parameter values as 
given in the caption of Fig. [S] 



points are generated. The latter hardly occurs at lower 
concentrations (compare Figs, [s] and [9] at time t = 8000). 
We have to note that the infiuence of the underlying 
square grid becomes rather obvious at the later numerical 
times. 

In order to describe the process of coarse-graining to 
the longer threads more quantitatively, we analyzed the 
kinetics of the end points. We refer to these points as 
defects in the following. There are several qualitatively 
different ways defects can evolve. Defects can appear 
together with the threads when the latter form from the 
homogeneous disordered state. This process mainly takes 
place at the very early stage of our numerical calculations 
and will not be considered. Apart from that, defects can 
appear when crosslinking points disconnect (inversion of 
the process shown in Fig.[7|^a)-(c)). Two defects can unite 
to one defect when threads become very short and the 
end points of the same thread meet. We have checked 
that this process does not play a major role in what 
follows. Defects further vanish by forming crosslinking 
points (in analogy to Fig. [7]^a)-(c)), or by connection of 
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FIG. 10: (Color online). Number of defects (end points of the 
threads) Ndef in the stripe patterns as a function of numerical 
calculation time t in a semilogarithmic plot. As marked by 
the straight line, the defect number decays exponentially in 
a broad regime from t — 800 to t = 6400. The inset shows 
the absolute displacement distance of the defects Ar per time 
interval At = 20 and averaged over the whole numerical sam- 
ple. Dashed vertical lines mark the interval used for the line 
fit. Parameter values as given in the caption of Fig. [S] 



two threads to one at the end points. The latter process 
dominates particularly at lower mean concentrations ip. 

Fig. 10 shows the number of the defects N^ef on a 
calculation grid of 512 x 512 lattice sites as a function 
of time. The mean concentration was = 0, and the 



other parameter values were chosen as mentioned in the 
caption of Fig. [Sj Directly after the initial quench from 
the homogeneous state to a temperature given by = 5, 
the number of defects decreases rapidly. As given by the 
straight line in the semilogarithmic plot, we then found 
a broad regime from about t = 800 to t = 6400, where 
the defect number decreases exponentially. At later times 
the annihilation of defects seems to saturate. Apart from 
that, we tracked the motion of the defects. In the regime 
of exponential decay, the average distance of motion of 
the defects (Ar) was relatively constant with about 0.34 
lattice sites per time interval At = 20. The correspond- 
ing standard deviation was less than 0.04 lattice sites. 
We have plotted the time dependence of (Ar) in the in- 
set of Fig. [lOj There, the dashed vertical lines mark the 
interval that we used for the line fit in the main part of 
the plot. 

We performed three independent numerical runs with 
different initial spatial Gaussian fluctuations around the 
homogeneous state. Out of these three runs, we chose 
the one where the defect number saturated the latest. 
The other two samples with earlier defect number sat- 
uration did not reveal a well-defined intermediate expo- 
nential regime. 

At the end of this section we point out that variations 
of the parameters /3, and b allow for further, quali- 




FIG. 11: Time series of the kinetic evolution of Eqs. ([T]), 
([2]), and ([7|), with a mean concentration -0 = 0. Parameter 
values were set to -i^ = 8, d = 15, /3 = 0.8, and b — 10, 
the grid size was 200 x 200 lattice sites of distance c/x = 1, 
and time steps were of size dt = 0.0008. First an interlaced 
structure of thinner threads form. These become unstable and 
partially coarse-grain to threads of thickness d. Then parts 
of the threads grow thicker than d, and finally a bubble-like 
structure of homogeneous bubble size and wall thickness d 
emerges. 



tatively different structures. We give one example in the 
There, the mean concentration was 



timeseries of Fig. 11 
set to = 0, and the other parameters were chosen as 
I? = 8, d = 15, /3 = 0.8, and b = 10. For higher thick- 
nesses d, we usually observe that an interlaced structure 
of thinner stripes forms from the spatially homogeneous 
state (first picture of Fig. 11 ). These stripes then become 
unstable with respect to the formation of the threads of 
thickness d. In the case of the current set of parameters, 
the threads then partially thicken to threads of a cross 
section higher than d (third picture of Fig. 11). Finally, 
the latter transform into a foam-like texture with walls 
of thickness d and uniform bubble size. The macroscopic 
phase separation and coarse-graining into blobs as ob- 
served in Figs. |2j |3j and [6] does not occur for this set of 
parameters. 



V. POSSIBLE APPLICATIONS 

In this section we discuss two possible applications of 
our energy density ([T]) and ^ to model systems in the 
field of soft condensed matter. In the first case, ?/^(r) is 
interpreted as a density field, which is therefore of the 
same spirit as is the phase field crystal model [13 for 
periodic systems. In the second part, ?/;(r) should rather 
be interpreted as a concentration field. This can be seen 
in analogy to the concentration field description of block 
copolymer systems [1 in the periodic case. 

From the density field point of view, as already men- 
tioned in the introduction, it appears natural to use 
Eqs. ([1]) and ^ for a density field characterization of 
polymeric systems. tp{r) is then connected to the time 
averaged probability of finding a segment of a polymer 
chain at position r, in analogy to the phase field crystal 
approach [19 . Our model then describes polymeric sys- 
tems on the length scale of a segment and on diffusive 
time scales. 

Two steps on this way to a realistic microscopic con- 
tinuum characterization of polymer systems seem to be 
important. First, we need to work out a more isotropic 



10 



discretization scheme than the one presented in the ap- 
pendix. This is necessary to decouple the orientation of 
the threads from the orientation of the calculation grid. 
Second, several topological properties of polymer systems 
can only arise in three dimensional space. This concerns 
topological entanglements that are not possible in two 
spatial dimensions. Special, yet illustrative, examples in- 
clude mixtures of ring and linear polymers, where the 
rings can be threaded by the linear chains [28l [29] , and 
slide-ring gels ^3Qrr32j . It will therefore be interesting to 
generalize our approach to three dimensional space. 

These points are important, e.g., for rheological consid- 
erations. Here, we leave them for future work. Instead, 
we now concentrate on the kinetics of the reaction pro- 
cess illustrated for example in Fig. [9] We note that the 
process of phase separation and formation of the stabi- 
lized stripes reflects several characteristics of step-growth 
polymerization reactions [33, 34 . In these reactions, all 
reactants/monomers can participate in the process from 
the beginning and are quickly incorporated into short 
chains. The chains can grow at both ends, and short 
chains can connect to longer ones. This means that the 
chain length continuously grows with increasing reaction 
time. Long reaction times are necessary to obtain high 
degrees of polymerization (long chains). All of these fea- 
tures are reflected by the process referred to in Fig. [9j 

As common in the literature, we now denote by x the 
number of units/monomers that a chain consists of. A^o 
gives the total number of units/monomers in the sample, 
whereas N denotes the total number of molecules of all 
sizes at a specific reaction time. Then the expression 
p = (A'o — N)/No describes the extent of reaction at that 
specific reaction time. As shown by Flory, the number 
Nx of molecules that consist of x units/monomers is given 

by 



TV, = iV(l-p)p^ 



(8) 



We now compare this distribution to the results from 
our numerical solution of Eqs. ([l]), ([2|, and ([7|) for = 
—0.4 (see Fig. [9|. For this purpose, we consider each 
pixel of > as holding a unit. To determine N^^ we 
extracted at each time the number of threads that consist 
of X connected pixels of > 0. N is given by the total 
number of threads (connected objects of > 0) at this 
time. In order to identify the total number of all units 
A^o, we counted the total number of pixels of > at a 
late reaction time, which was at t = 800. 

The distributions were determined from nine separate 
numerical samples of different initial fluctuations around 
the spatially homogeneous state. Each sample was of size 
512 X 512 grid points, the parameter values were chosen 
as given by the caption of Fig.|9] For small values of x, we 
expect a significant deviation of from our numerical 
results. The reason is that our process of phase separa- 
tion does not support objects of a diameter smaller than 
d. Also, diffusion can be important for small objects, be- 
fore the process of chain connections becomes dominant. 
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FIG. 12: (Color online). Time evolution of the chain/thread 
length distribution Nx^ where x gives the number of units in 
the chain/thread. The numerical data points were obtained 
from nine independent numerical samples at the times indi- 
cated on the upper right. Corresponding parameter values are 
given by the caption of Fig. [9] Each point is an average over 
20 sizes. The lines follow from the Flory chain length distri- 
bution for step-growth polymerizations as given by Eq. (jsj. 
The necessary parameters were extracted from the numerical 
results. 



Fig. 12 compares the results from our numerical solu- 
tions to the predictions of Eq. ([8|. The numerical data 
points are averaged over 20 object sizes. We chose this 
size of averaging because it corresponds to the minimal 
supported object size of 7r((i/2)^. Clearly, we can see the 
deviations outlined above. Considering, however, the fact 
that there is no adjustable parameter, the description by 
Eq. ([8| works reasonably well. It allows an approximate 
connection to step-growth polymerization processes. 

In the second part of this section, we explain that a 
concentration field interpretation of Eqs. ([T]), ([2|, and 
^ can be seen as a phenomenological approach to qual- 
itatively describe micellar, vesicular, and bilayer mem- 
brane systems. For this purpose, we think, e.g., of sur- 
factant molecules in an aqueous solution. The surfactant 
molecules consist of a hydrophilic head group and a hy- 
drophobic tail. We assume that the mass density of the 
hydrophilic and hydrophobic parts is identical. Then we 
can characterize the pure hydrophilic components (water 
and hydrophilic head groups) by a concentration = — 1, 
and the pure hydrophobic parts by t/^ = -hi. 

The phase separation for high enough values of cor- 
responds to the separation of hydrophilic and hydropho- 
bic components. For example, in Fig. [ija), the spheres 
would correspond to micelles, and the stripes to bilayer 
membranes. The hydrophilic head groups are on average 
located on the dark side {^Ij ^ —1) oi the dark-to-bright 
boundaries. On the contrary, the hydrophobic tails point 
into the center of the bright regions (?/^ ~ +1). 

Since hydrophilic and hydrophobic parts are connected 
to each other, the width of the hydrophobic region can at 
maximum be twice the average length of the hydrophobic 




tails. We therefore roughly identify the thickness d of 
the bright regions in Fig. [TJa) as twice the length of the 
hydrophobic tails. The /3-term in Eq. ^ determines the 
degree of depletion of hydrophobic groups around the 
regions of hydrophilic head groups. 

It may be possible to relate processes similar to the 
ones shown in Fig. [7] to the formation of a vesicle. In the 
early stages of the numerical calculations, or throughout 
the calculations of section \P7\ the threads would have 
to be interpreted as wormlike micelles. The higher cur- 
vature energy density at the ends of wormlike micelles 
can be mimicked by our energy functional, as, e.g., de- 
picted in the upper picture of Fig. [4] In the following, 
however, we will restrict ourselves to qualitative hydro- 
dynamic considerations. More precisely, we investigate 
the flow behavior of vesicular structures in simple shear 
and pipe flows. 

For this purpose, we coupled Eqs. ([T]), ([2|, and ^ with 
the Navier-Stokes equation, which leads to ^ 



'dt 

'di 



-V • VV^ + LA 



(9) 



-Vp - V • Vv + T^Av - - V^)V^. (10) 

dip 



Here, we have assumed incompressibility. This makes the 
continuity equation reduce to V • v = 0. L is the same 
parameter as in Eq. and rj denotes the viscosity. The 
(constant) mass density p has been scaled out. 

We set the parameter values to = 5^ a = 3^ d = 
15.3, f3 = 0.3, b = 10.2, and r] = 100. In the case of 
simple shear flow, we set L = 10, whereas we chose L = 
20 in the case of pipe flow. The numerical calculation 
grid was of size 256 x 128 sites with a lattice distance of 
dx = 1, and the size of the time steps was dt = 0.00005. 
We chose periodic boundary conditions at the left and 
right edges of our rectangular calculation grid. On the 
contrary, noflux boundary conditions were imposed at 
the top and bottom boundaries. The flow direction will 
thus be horizontal. 

To enforce steady simple shear flow, we set the hori- 
zontal velocity equal to —2 at the upper boundary. In all 
other cases, both velocity components were set equal to 
zero at the top and bottom boundaries. To induce steady 
horizontal pipe flow, we added a constant gravity- like vol- 
ume force g = (0,^) to Eq. (10), where ^ < in our case. 
The horizontal flow will therefore be directed towards the 
left edge of the calculation grid in both geometries. We 
used a MAC scheme on a staggered grid to numerically 
integrate Eqs. ([9| and (10) forward in time [37 . The 
intermediate Laplace equation resulting for the pressure 
from the incompressibility condition V • v = was solved 
through a spectral method [38^ . Third and first order up- 
wind schemes were used to discretize the advective term 
in Eqs. (9| and (10), respectively. 

In order to start our investigations, we provided an 
annulus of = and thickness d in an environment of 
ip = —1. We let it relax to a spherical vesicle-like object 
by numerically iterating Eq. ([2| for a time interval At = 




FIG. 13: (Color online). Simple shear flow of a vesicle char- 
acterized by Eqs. ([9| and (10). An approximately spherical 
vesicle was provided as an initial condition. At time t — 5 the 
shear was started by setting the horizontal velocity compo- 
nent at the upper boundary to a nonzero value, pointing to 
the left, the effect of which is clearly visible at time t — 27.5. 
The bottom row shows the vesicle after having crossed the lat- 
eral periodic boundaries once (time t = 200) and twice (time 
t = 455). Pink color: ip > 0.8, red color: ip > 1. Parameter 
values were set to = 5, a = 3, d = 15.3, /3 — 0.3, b — 10.2, 
7] = 100, and L = 10. The grid size was 256 x 128 sites with a 
lattice distance of = 1, and the size of the time steps was 
dt = 0.00005. 



40. These vesicles were then used as an initial condition 
for the numerical investigations of Eqs. (I9|) and (10). The 



objects that we describe by this sort of concentration field 
approach have to be classified as fluid vesicles. That 
is, viscous flow of the constituents is possible within the 
vesicle membranes. No elastic shear forces occur within 
the local membrane planes, which would be the case for 
elastic membranes. 

We give an example for the obtained behavior of a 
vesicle under shear flow in Fig. 13 The first picture shows 
the vesicle at the time when the shear flow is started by 
setting the horizontal velocity component to a nonzero 
value at the upper boundary. In the second picture, the 
effect of this shear deformation starts to become visible. 



The bottom row of Fig. [13] depicts the vesicle as it is 
strongly deformed by the flow field of steady shear. We 
show its state after having crossed the periodic left-right 
boundary once and twice. 

The long axis of the stretched vesicle approaches a 
steady orientation with respect to the flow direction. 
This behavior is well known from numerical studies, 
e.g. solvent multiparticle collision dynamics combined 
with coarse-grained membrane models [39l |40], as well 
as experimental observations [42]. It is referred to 
as the tank-treading regime, because it is usually accom- 
panied by a rotational motion of the vesicle membrane 
around its interior. 

We must remark in this context that our characteri- 
zation in its current implementation does not reproduce 
an essential feature of the experimental systems. Real 
vesicles are practically impermeable for the surrounding 
solvent on experimental time scales. Generally, in field 
descriptions, this can be achieved by setting the mem- 
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brane velocity equal to the local fluid flow velocity [43]. 
In our case, this corresponds to letting L ^ in Eq. ([9|. 
Then the concentration fleld describing the location of 
the membrane is advected by the flow velocity only and 
does not diffuse. However, this limit is not accessible by 
our current numerical scheme and is therefore beyond the 
scope of this work. 

Analyzing the flow fleld around the vesicle, we do, how- 
ever, observe the tendency of tank-treading along the 
vesicle membrane. Fig. p!4|^a) shows the velocity proflle 
when the membrane is exposed to a steady shear flow. 
The situation corresponds to the last picture of Fig. [13] 
We flnd the flne structure of the flow velocity fleld by sub- 
tracting the linear shear proflle, as depicted in Fig.[T4]^b). 
This reveals how the membrane influences the velocity 
proflle of the surrounding fluid. Two vortices form at 
the front (left) and two at the back (right) of the vesi- 
cle (the bottom one at the front is not as clearly visible 
as the others). Within these two groups of vortices, one 
vortex adds to aligning the elongated vesicle along the 
flow direction, the other to aligning it perpendicular to 
the flow fleld. This competition is reflected by the flnite 
inclination angle of the stretched vesicle. Furthermore, 
the membrane connects the flow in the upper half of the 
channel to the flow in the lower half and slightly evens 
the different flow velocities. In Fig, p^b), this becomes 
clear from the velocity fleld predominantly pointing to 
the right in the upper part, whereas it is oriented toward 
the left in the lower part. Generally, such a balancing 
effect between the areas of different flow velocities can 
be connected to the tank-treading motion of the mem- 
brane. Fig. 14 ^c) displays the flow fleld of the membrane 
in its rest frame. A rotational tank-treading-like motion 
of the membrane is revealed by rescaling the horizontal 
component of velocity by a factor of 0.01 with respect to 
the vertical component. 

Apart from the tank-treading motion, at least two 
other regimes were identifled numerically and experimen- 
tally [39H42] |44] . If the inclination angle of the long axis 
of the elongated vesicle is not steady but varies in time, 
the motion is called tumbling. In addition, a more ir- 
regular dynamic mode of motion was identifled: it in- 
cludes tumbling, features dynamic shape deformations 
away from the purely elongated state, and was named 
swinging mode. 

Mainly two parameters were found to determine the 
mode of motion. First, this is the reduced volume. It 
measures the deviation from a spherical shape by relating 
the surface to the volume of the vesicle. The second 
parameter is the ratio between the fluid viscosities inside 
and outside the vesicle. Imposing impermeability of the 
membrane is crucial to control these parameters. It will 
therefore be an important step for future work on our 
description. 

We close this section by pointing out our qualitative 



results for the pipe flow geometry. Fig. [15] shows an ini- 
tially spherical vesicle that has a size of the order of the 
channel diameter. When the flow starts, the vesicle is de- 
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FIG. 14: (Color online). Flow field v corresponding to the 
case of simple shear in Fig. |13| at time t — 455. The red color 
indicates the flow field within the membrane (-0 > 0.5), the 
blue color the flow field within the remaining solvent, (a) 
Steady shear fiow acting on the membrane, (b) Subtraction 
of the linear shear profile reveals the fine structure and dis- 
tortion of the fiow field due to the presence of the membrane. 
Four vortices form, two at the front and two at the back of 
the membrane (the bottom one at the front is not as clearly 
visible as the other ones). Competition within the group of 
the two front vortices and within the group of the two vortices 
at the back is refiected by the finite inclination angle of the 
membrane, (c) Tendency of tank-treading in the rest frame of 
the membrane. For clarity, the horizontal velocity component 
is rescaled by a factor 0.01 w.r.t. the vertical one. 
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FIG. 15: (Color online). Time series of a vesicle exposed 
to pipe flow. First, the vesicle folds into a parachute shape, 
following the Poiseuille flow profile. Then it tends to elongate 
along the flow direction. Parameter values and color code 
as given by the caption of Fig. [13] except for L = 20 and 
9 = -0.4. 
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FIG. 16: (Color online). Two vesicles in a pipe geometry. As 
time proceeds, the rear vesicle grows on cost of the front vesi- 
cle. Parameter values and color code as given by the caption 
of Fig. [13] except for L = 20 and g = —0.5. 



formed according to the profile of the Poiseuille-like flow. 
It folds into a parachute-like shape. At longer times, we 
find a transition to a shape that is elongated in the flow 
direction (compare Ref. p5^). 

Finally, we put two vesicles into the channel, as de- 
picted by Fig. [T6| Initially, they show the same behavior 
as the single vesicle in Fig. [15] However, as time pro- 
ceeds, the front vesicle looses material that is collected 
by the rear one. Consequently, the first vesicle shrinks, 
whereas the second one grows. It would be interesting to 
know whether a corresponding process can be observed 
in an experimental system. 



nally we observe a macroscopic phase separation. The 
scaling behavior of the latter coarse-graining process is 
similar to the case of the Cahn-Hilliard model. We could 
omit the macroscopic phase separation by adding expres- 
sion ([7]) to the energy density. It stabilizes the localized, 
non-periodic thread-like stripe structures. In this case, 
we observe a connecting process of different threads to 
reduce the number of high-energetic end points of the 
threads. Our analysis revealed a regime during which 
the number of end points of the threads decays exponen- 
tially. Finally, we discussed two possible applications of 
our model in the field of soft matter. On the one hand, 
this was a density field description of polymers, where 
for the moment we concentrated on the chain length dis- 
tribution function during the formation process. On the 
other hand, we focussed on the hydrodynamic proper- 
ties of vesicular structures. Simple shear and pipe flow 
geometries were qualitatively addressed, where our de- 
scription reproduces generic features of the flow behav- 
ior. 

As pointed out, several aspects remain for future in- 
vestigations. We have indicated that the variation of the 
various parameter values involved in the description leads 
to further qualitatively different textures. A complete 
phase diagram needs to be established. Next, to allow 
more quantitative studies, our numerical discretization 
scheme should be improved. For example, a further de- 
coupling from the underlying grid directions might be 
achieved through an implementation on a hexagonal lat- 
tice. Furthermore, a generalization to three spatial di- 
mensions would be desirable. Lastly, we have followed 
a phenomenological approach so far. It would be inter- 
esting to investigate whether our functional can be re- 
lated to more microscopic approaches through appropri- 
ate coarse-graining, and which changes to the functional 
form this would imply. 
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VI. CONCLUSIONS 



Appendix: Discretized numerical evaluation of the 
functional derivative and kinetic equation 



We have introduced a simple nonlocal energy density 
in Eq. ([T]) that was discretized for numerical investiga- 
tions as explained in the appendix. Starting from the 
disordered (spatially homogeneous) state, this potential 
describes the formation of non-periodic localized thread- 
like structures. These structures are metastable, and fi- 



In this appendix we describe our numerical evaluation 
of the functional derivative in Eq. ([2|. The form of the 
corresponding energy density is given by Eq. ([T]). 

Due to the nonlocal contribution in Eq. ([T]), it is not 
straightforward to perform the functional derivative ana- 
lytically. We use the following definition of the functional 
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derivative: 



1 r 



= lim 



^WrO + 6J(r'-r)}-^{^(rO} 

(A.l) 

This is the correct starting point for the transformation 
from the analytic continuum picture to the discrete nu- 
merical calculation. 

The numerical lattice is of finite size and defines a dis- 
crete system. We use a rectangular domain on a square 
lattice of x Ny lattice sites with periodic boundary 
conditions. The lattice sites are indexed by pairs (i, j), 
where z = 0, . . . , A^^^ — 1 and j = 0, . . . , A^^ — 1. Conse- 
quently, the continuous functional form J^{?/^(r)} trans- 
forms to a simple multivariate function T [^l^ij] of x Ny 
variables ipij. 

In analogy, the Dirac delta function 5{v' — r) is now 
given by a product of Kronecker deltas Sudjm^ where (i, j) 
and (/,m) correspond to the positions r and r', respec- 
tively. Therefore, the functional derivative in Eq. (A.l) 
transforms to a simple partial derivative. 



lim - 



]-T[^lm]l (A.2) 



This is the expression that we use to numerically evaluate 
the functional derivatives. 

Furthermore, we discretize the gradient field in Eq. ([T]) 

as 



dx 



dx 



(A.3) 



We can motivate this choice by considering a quadratic 
gradient contribution |(VV^ |r)^ in the energy density, as 
it occurs, e.g., in the energy functional for the Cahn- 
Hilliard model On the one hand, analytically, 

the functional derivative of this term follows immedi- 
ately as — A?/^ |r. On the other hand, the discretiza- 
tion of the term ^iVijj |r)^ according to Eq. (A.3) reads 

/dx^. 



\ - "^ijY + - V^ijO^J /dx^- We ca n per - 

form the functional derivative as described by Eq. (A.2). 
Equating both results at position | 



At/; \ij = 



yields 



dx'^ 



(A.4) 

As we can see, we correctly obtain the discretization of 
the Laplace operator as given by the stencil 






1 





1 


-4 


1 





1 






(A.5) 



This discretization of the Laplacian introduces 
anisotropy [46]. Mo re sophisticated discretizations 
than expression (A.3) must be worked out for the 
gradient field when a fully isotropic discretization of 
the Laplacian operator should be obtained from the 



functional derivative. This is beyond the scope of the 
considerations in this paper. 

In order to write our results in a compact way, we 
define an abbreviation for the discretized nonlocal sub- 
scripts I I ^ vy; I at site (i, j): 



\xpij,ypij 



l-(r+<i^)|,.y.(r+<i^)|^^.- 



(A.6) 



The pair {xp., 
script Ir+rf 



liv^ll 



yPij) indexes the lattice site that the sub- 
points to (the "p" stands for "pointer"). 



Here, the gradients in the subscript are discretized ac- 
cording to Eq. (A.3). 



On the one hand, the value tpij contributes to the value 
of the energy density at the grid site (i,^). This comes 
from the local "d-teim and from the local term \/tp\r in 
the square brackets of Eq. ([T]). It is crucial to note that, 
on the other hand, tljij contributes to the value of the 
energy density at all other grid sites (/,m) that point to 
{ij) through the subscript Upim^ypim = Ir'+d^- The 
latter part comes from the nonlocal term \/th\ , ^ vy. in 

the square brackets of Eq. ([T]). 

Consequently, at each time step, we must determine 
for each lattice site (i, j) all other lattice sites (/,m) that 
point to (i, j); i.e. all other lattice sites (/,m) satisfying 
the condition {xpim^ypim) = {hj)- In the supplemental 
material [26 we present a scheme to perform this task at 
an overall computational cost of linear order, 0{NxNy). 
We used dynamic list structures for this purpose. 

Following the procedure outlined in Eq. (A.2), we ob- 



tain the contributions to the functional derivative at site 
(z, j). From the i^-term in Eq. ([T]), we directly find 



(A.7) 



Next, at site (^, j), the local term VV^|r leads to a con- 
tribution 



'^Xpim-\-'^,yPlm '^Xpim.ypi 



l,m 

x(^Z+l,i(^mj — ^li^mj) 

H~ (V^^jm+l "^Im H~ "^^Xpim^yPlm-^-^ "^Xpim^yPlrr 

X{5ii5m+lj - 5ii5mj) j. (A. 

For the nonlocal term V'^lr \ d ' obtain a con- 



"liv-^ii 



2a 

dx'^ 



Lm 



tribution 

^ ^ I {^l-\-l,m ~ ^Im + "^Xpim+l.yPlm ~ ^Xpim,yPlm) 
^{^Xpim-\-l,i^yPlm,j ~ ^Xpim,i^yPlm,j) 
+ {^l,m-\-l ~ ^Im + ^Xpim,yPlm-\-l ~ ^Xpim,yPlm) 
^{^Xpim,i^yPlm+^j ~ ^Xpim.i^yPlmj) | (^'9) 
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to the functional derivative at site (i,^). 

Finally, we have to add the contribution arising from 
the term with the coefficient /3 in Eq. It turns out 
that the gradients in the subscript |^ ^ must be dis- 



liv^ll 



cretized differently than in Eqs. (A.3) and (A.6). Oth- 
erwise the discrete nature of the numerical calculation 
leads to a systematic artificial drift. 
In analogy to Eq. (A.6), we define 



\xp(3ij,yp(3i. 



livv-ll 



(A.IO) 



Here, however, the gradients in the subscript are dis- 
cretized according to 



2dx 



2dx 



(A.ll) 



instead of Eq. ( |A.3[ ). For each site (/,m) that satisfies 
the condition {xp/3im,ypf^im) = {hj) we obtain a contri- 
bution 



2p{^,j + 1) (A.12) 



to the functional derivative at site (^, j). 

We present a schematic implementation of the above 
results using a C-like pseudocode in the associated sup- 
plemental material ^26] . 
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